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ABSTRACT 

Professor  R.  D.  Rlchtmyer  has  described  a  finite  difference 
method  for  the  computation  of  hydrodynamlc  flows  which  contain 
a  shock  [3j^]-   This  method  uses  the  Bulerlan  form  of  the  hydro- 
dynamic  equations,  is  explicit.  Is  of  second  order  accuracy, 
and  Is  based  on  shock  fitting  rather  than  the  introduction  of 
artificial  viscosity.   This  paper  describes  the  result  of 
numerical  computations  using  a  modification  of  this  finite 
difference  method.   The  method  is  applied  to  a  one-dimensional 
problem  for  which  a  solution  can  be  computed  by  solving  an 
ordinary  differential  equation.   Therefore  we  are  able  to  de- 
termine the  accuracy  of  the  method  for  this  problem. 
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NYO-9605 

NUMERICAL  COMPUTATION  OF  HYDRODYNAMIC 
PLOWS  WHICH  CONTAIN  A  SHOCK 

1  .   Introduction . 

Our  objective  is  to  test  the  method  of  Rlchtmyer  for 
stability  and  accuracy  by  using  It  to  compute  the  one- 
dlmenslonal  flow  of  a  shock  down  a  tube.   The  shock  fitting 
scheme  is  a  slight  modification  of  that  proposed  by 
Rlchtmyer  [5,4]  .   Two  finite  difference  schemes  for  the  flow 
behind  the  shock  are  described  and  tested.   We  will  describe 
the  problem  used  for  the  test  calculation,  the  finite  difference 
methods,  and  then  the  results  of  the  calculations  performed  on 
the  IBM  7090. 


2 .   Description  of  the  Problem. 

The  problem  for  which  we  will  compute  the  solution  is  that 

of  a  shock  moving  down  a  tube  with  a  variable  density  cold  gas 

ahead  of  the  shock.   This  problem  was  suggested  by 

Professor  Peter  Lax.   It  has  the  advantage  of  providing  an 

exact  solution  with  a  shock  of  variable  strength.   The  problem 

and  its  solution  is  due  to  J.  Keller  [l].   We  denote  the  distance 

along  the  tube  by  x   and  the  pressure,  density  and  velocity 

ahead  of  the  shock  by   P  ,  p  ,   and  u  ,   respectively.   Keller 

0*^0         o 


showed  that  if   P  =  u 

o    o 


p   =  2F^(7+l)~^R^x  ^ 
'  o     o        o 


0  and 
2a 


--^  r2  +  a^Rj'^^   X 
I-7   o     2  o 


H_       2(7-1) 
7+1 


3 

2 
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then  the  flow  behind  the  shock  can  be  determined  by  solving 
the  following  differential  equation 


j"(t) 


2«^    (1-r) 


+  a^ 


1 
2 


In  the  above  equations  ^   ,    B.   ,    a-,  ,      and   ap  are  arbitrary 
constants  and  7   is  the  ratio  of  specific  heats.   If  we 
define  the  function  F(x)  by 

,2 


F(x)  =  F  X 

o 


2«lRo  ^    2 
1-7 


1_ 
'2 


then  the  flow  behind  the  shock  is  given  by 

p(t,x)  =  -  ^'(^y^i^^^  p(t,x)  =  j-^(t)F(xr^(t)) 


a-jX 


-1 


u(t,x)  =  xj^(t)j  ^(t) 


j(t) 


7+1 
2 


In 


The  shock  position  is  given  by   R(t)  =  R 
our  computations  we  use  the  values   P  =  ^0,    R  =  O.O868, 
a-,  =  -  1,  ttp  =  0.195^   and  7  =  3-   In  figure  2  graphs  of 
p(0,x),  p(O^x),   and  u(0,x)   are  shown.   The  function   j(t) 
is  evaluated  by  a  Runge-Kutta  routine;  or  If  7  =  J> ,    then  the 
differential  equation  for   j(t)   Is  solved  by  an  explicit 
formula.   The  flow  Is  computed  by  the  finite  difference  scheme 
between  x  =  a   and  the  shock.   The  boundary  conditions  at 
X  =  a  are  obtained  from  the  exact  solution  described  above. 
The  flow  is  then  determined  from  the  usual  hydrodynamlc 
equations  behind  the  shock  and  the  Hugoniot  relations  for  the 
shock. 
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3-   The  Finite  Difference  Equations . 

We  denote  the  energy  and  momentum  by  e   and  m, 
respectively  and  make  the  definitions 


V  = 


P 
e 

m 


f(v)  = 


m 


7me    7-1  m 


3 


2   P 


The  flov;  equations  in  conservation  form  are 


^v  ,  ^f    rs  ^v  ,  .  ^v     „ 

—  +  —  =0   or  —  +A  —  =0 


^t   dx 


^t 


bx 


where   A   is  the  Jacobian  of   f (v) .   Two  finite  difference 
schemes  are  tested  on  this  problem.   The  first  is  the  scheme 
of  Lax  and  Wendroff  [2]  and  the  second  is  an  Iterative  scheme 
to  be  described  below. 

The  difference  equations  used  by  Lax  and  Wendroff  are 
the  following: 


n+1    n 

V.     =  V. 


,n 


-(f"    -  f"   )  +  - 
2^1+1     1-1''  ^  2 


2  r 


,n  /  „n 


,n ' 


,n 


,n 


^1^  1+1     1^    ^1-1^  1     1-1 


n 


where   A  =  7—  ,  v.  =  v(t  ,x.)   and   similarly  for   A.   and   f; 
Ax  '   1      n   1  ^  1        1 

These  equations  are  based  on  an  expansion  of   v(t,x)   in  powers 
of  At  about  [t      +   -p^J  •   Using  the  flow  equations,  we  find  the 

first  two  terms  of  the  expansion  to  be 

,1 


av 


'2        /-..Nn        (^^(a|^) 


n 


A  stability  analysis  by  von  Neumann's  method  Indicates  that 
this  scheme  Is  stable  provided  that   A  <   -i — t-  where   a   Is  the 
largest  eigenvalue  (in  absolute  value)  of  any  of  the  matrices 
A_.  .   This  analysis  Is  based  on  the  linearized  equations  and 

Ignores  the  Influence  of  the  boundary  conditions.   The 

n  "^ 

truncation  error  In  v.   Is   O(Ax^). 

We  denote  the  mesh  point  Immediately  behind  the  shock  by 

X.    and  denote  the  shock  position  by  R  =  R(t  )   (strictly 
IS  ri 

speaking   Is  =  ls(n),   see  figure  l).   The  above  difference 

equations  can  not  be  used  at  x.    since  x.  , -,   does  not  lie 
^  IS  is+1 

behind  the  shock.   In  order  to  compute   ^f/^x  with  second 
order  accuracy  we  must  use  the  values  of   v  at  three  points. 

We  can  use  the  points   x.   o?  ^-   ij   ai^d  x.    or  x.^  -,, 

n 
X.    and  R  .   The  results  of  the  computation  show  that  the 

first  choice  is  somewhat  better,  although  the  difference  is 

not  significant.   The  difference  equations  at  x.    are, 

using  the  first  choice, 

v'^+1  =  v'^   -  -(f     -  ^f     +  3f   ) 
is    ^is    2^  is-2   ^^is-1  ^  -^  is^ 


a2  r 


+  2 


n  ,  n   _   n    )  _  a^    (^n    _   n    ^ 
^Is^  IS     is-1^    ^is-l^^ls-1     is-2^ 


The  second  finite  difference  scheme  is  an  iterative  scheme 


defined  by 


0  n 

CO.  =  V. 

1  1 

n+1  p 

V.  =  tQ^ 

1  1 

,s+l  n    A 

CO.  =  V.   -  TT 

1  1    4 


fK+i'  -^K-i'  +f<-?+i'  -f'-;-i> 
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This  scheme  Is  described  In  more  detail  elsewhere  [5]-   The 
scheme  Is  stable  If  A  <  2/|a|   and  p  =  3,4,7,8,...    ana 
unstable  otherwise.   (a  Is  the  largest  eigenvalue, ' in  absolute 
value,  of  any  of  the  matrices  A..)   When  1  =  is  we  must 
modify  the  difference  quotients  in  the  same  manner  that  they 
were  modified  in  the  Lax-Wendroff  scheme.   If  p  ^  2,   then 
the  truncation  error  is   0(Ax''). 

4.   The  Shock  Fitting  Method. 

The  shock  fitting  method  is  essentially  that  proposed  by 

R.  D.  Rlchtmyer  [3].   We  let  p  ,e  ,  m   denote  values  Immedi- 

s   s   s 

ately  behind  the  shock  and  let  V  denote  the  velocity  of  the 
shock.   Thus  there  are  four  unknowns  at  the  shock.   The 
Hugoniot  relations  provide  three  equations  relating  these 
unknowns  at  the  shock.   The  three  flow  equations  must  provide 
the  additional  relation.   Only  one  of  the  three  characteristics 
at  the  shock  lies  behind  the  shock.   Rlchtmyer  proposed  the  use 
of  the  corresponding  characteristic  equation  as  the  fourth  re- 
lation.  He  also  conjectured  that  the  use  of  a  characteristic 
equation  corresponding  to  a  characteristic  lying  ahead  of  the 
shock  would  lead  to  instability.   This  conjecture  is  tested 
and  it  is  correct. 

This  characteristic  equation  can  be  written  as  follows 

hv  bf  .bv  dv 

where   A  is  the  Jacobian  of   f(v)   and   r   is  the  eigenvector 
of  A  corresponding  to  the  eigenvalue   a  =  u  +  c  (c  =  s/yp/p) • 
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The  Hugoniot  relations  can  be  written  as 


V  (v   -  V  )  =  f (v  )  -  f (v  )  =  f   -  f  , 
s    o       s       o     s    o 


The  shock  fitting  is  done  by  iteration.   It  is  the  same  regardless 

of  which  scheme  is  used  for  the  flow  behind  the  shock  (that  is, 

Lax-Wendroff  or  iteration)  .   We  ass^ome  that  approximations  for 

Vg    and  V     are  known,  then  we  compute  corrected  values 

v?"^^  =  v^"^^  +  Av   and  V^"^^  =  v"^"^^  +  AV  as  follows, 
s      s 

If  we  denote  the  directional  derivative  along  the  shock  by 

D 


/Dt ,   then 


Dv  _  ^_v   Tj  hv_ 


Prom  the  characteristic  equation  we  have 

Dv         /  „  -r-r  ^^v 

The  above  equation  in  finite  difference  form  and  the  Hugoniot 

Cb  Y^ 

relations  are    (  <  -^  >  must  be  evaluated  by  a  difference 

quotient,  which  will  be  defined  below) 

A  ,1.1      , 1    /   K  n-4- 

2./— n+1     n\      .,     2  /.   2    ,r  2-r\  >  ov  ) 

r„   '(v     -v)=-Atr    "(A    -V   I)  ? 3— f 

s     s      s  s     s      s     C  °^  J 

T7n+1/— n+1      \    „/— n+l\    „/   \ 
V    (Vg    -  v^)  =  f(Vg   )  -  f(v^). 

Expanding  these  equations  out  to  the  first  order  In  Av  and 

y-<  I  J- 

AV  we  have   (let   r  =  r^  ^,  J  =  Ta  -  Vl)^'^"'- 
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i'-€Mil 


1 

n    ^,    C  ^r)'"2   Cdv)'^^      n+1 

r.Av  =  r-v^  -  At  ^a  -  Vj^   ^  fe)^    "  ^'^s 


■'(' 


s    o  oj  Is    or 


Av  =  J~^\v{v      -  V  )  -  (f  -  f  )f'^'^^   +  AVJ"^(v^  -  v^ 


These  equations  are  easily  solved  for  AV  and  Av .   The 

\  ^^'^  ^  J-    •   J  4^      n+1   n  -rrn+l 

average  values  i     (  are  determined  from  v    ,  v  ,  V 


n 
and  V   and  are  therefore  Independent  of  Av  and  AV, 


^1         ^_^ 
The  value  of  ^  ^j^     is  ^^j^    =  ^ 


n+1    ,   .n 


The  latter  derivatives  must  be  computed  to  second  order  accu- 
racy in  Ax.   Two  methods  of  computation  are  tested.   The 
first  is  to  pass  a  second  degree  polynomial  through  the  values 
of   v"^"^    (or  v'"^)   at  the  points   x.   -,,  x.    and  r'^    (or  R*^). 

■^  lS-1     IS 

The  second  is  to  fit  a  second  degree  polynomial  by  least  squares 

through  the  values  of   v'^^   (or  v'^)   at  the  points   x.   „, 

x.   -,  ,  X.    and  R*^"*"   (or  r'^  )  .   The  first  method  is  satisfactory, 

lS-1     IS 

the  second  is  not  (see  the  discussion  of  the  results  below). 

Thus  we  have  an  iterative  method  for  computing  the  solution 
at  the  shock.   If  the  iterative  scheme  is  used  for  the  flow 
behind  the  shock,  then  one  iteration  at  the  shock  is  performed 
for  each  iteration  behind  the  shock.   If  the  Lax-Wendroff  scheme 
is  used  behind  the  shock,  then  the  solution  behind  the  shock  is 
computed  before  the  iterations  are  performed  at  the  shock.   In 
the  Lax-Wendroff  scheme  three  iterations  are  generally  used  in 
the  shock  fitting. 

We  must  still  describe  how  this  scheme  computes  the  solution 
at  a  mesh  point  which  is  overtaken  by  the  shock,  for  example  the 
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point   (t  ,-,,  X.  ,,)   In  figure  1.   It  would  obviously  be 
^  n+1    is+1 

awkward  to  use  the  flow  equations  to  compute  the  solution  at 
this  point;   therefore  Interpolation  Is  used.   A  second  degree 


at  the  points 


n+1 
polynomial  Is  passed  through  the  values  of   v 

X.      -,  ,    X.        and  R     and  this  polynomial  Is  used  to  Interpo- 

n+1 
late  for  v.   -,  .   This  Is  the  same  polynomial  used  to  compute 

/^^\n+l    ^^"^^ 

7  Sx")     ■   ^^  could  also  use  a  polynomial  obtained  by  a  least 

squares  fit.   However,  this  least  squares  fit  Is  also  unsatis- 
factory. 


The  value  of  At   is  changed  at  each  time-step  according  to 


the  formula  At  = 


pAx 


where   a  =  |u   +  c  I  ,  c^  =  \/y?   /p 


m 


m 


and  0  <    p  <  1.   Usually  we  take   p  =  0.9-   The  program  will 

stop  if  the  shock  overtakes  more  than  one  mesh  point  per  time-step 

5-   The  Results  of  the  Calculation. 

The  tables  below  show  the  maximum  percentage  error  in  the 
pressure  after   20,  40,   and  100   time-steps  for  various  values 
of  Ax.   Unless  otherwise  stated  the  least  squares  fit  was  not 
used  to  compute  the  derivatives  or  to  interpolate  for  the  over- 
taken mesh  point.   Also,  unless  otherwise  stated  At  =  0.9Ax/a  • 

Table  I.   Percentage  error  for  the  Lax-Wendroff  method. 

Ax     Time-steps 


20 


40 


100 


.05 

0.17 

0.52 

0.69 

.025 

0.022 

0.043 

0.098 

.0125 

0.0029 

0.0055 

0.013 
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Table  II.  Percentage  error  for  the  Iteration  method. 


Ax 

Tlme-st 

ep 

s 

No 

.  of 

20 

40 

100 

200 

Iterations 

0.05 

0.17 

0.32 

0.69 

-- 

3 

0.025 

0.025 

0.046 

0.098 

-- 

5 

0.012  5 

0,0049 

0 . 007  5 

0.014 

-- 

5 

0.05 

0.l6 

0.32 

0.69 

-- 

4 

0.05 

0.32 

0.41 

0.76 

1.20 

2 

Note  that  there  is  very  little  difference  in  the  accuracy  of 

the  two  methods.   Inspection  of  these  results  shows  that  the  error 

is  approximately   0(z\x  )   especially  for  the  Lax-Wendroff  method. 

■  Also  note  that  the  iteration  method  shows  no  sign  of  instability 

out  to  200  time-steps  when  only  two  iterations  are  used.   This  in 

spite  of  the  stability  analysis  which  indicates  Instability  for 

two  iterations.   This  is  probably  due  to  the  fact  that   At  =  0.9Ax/a^ 

where  the  maximum  allowable  value  of  At   for  the  iteration  method  is 

At  =  2zXx/a  .   If  the  larger  value  of   At   is  used  the  shock  will 
'    m  ^ 

overtake  two  mesh  points  per  time-step.   Other  experiments  with 
the  iteration  scheme  have  indicated  instability  when  At  =  1.9Ax/aj^ 
and  apparent  stability  with  At  =  0,9Ax/a^  when  only  two  iterations 
are  used  [ 5] • 


m 


Table  III.   Percentage  error  for  the  Lax-Wendroff  method  when 

the  difference  equation  for  v?    uses  the  values  of   v   at 

X.   -,,  X.  ,   and  R    (see  section  3)- 
is-1^   is^ 
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Ax     Time-steps 
20       kO 


0.05 
0 .  02  5 


100 


0.25 
0.031 


0.51 
0.043 


0.69 
0.098 


Comparison  with  table  1  shows  that  slightly  better  results 

are  obtained  using  the  values  at  x.   ^,   x.   -,  ,   and  x.    in 

^  is-2'    is-1         IS 

the  difference  equation  for   v.   .   In  this  case  the  values  of 

^  IS 

V   on  the  shock  are  never  used  in  the  computation  of  the  flow 
behind  the  shock.   When  the  shock  overtakes  a  mesh  point  the  values 
of   v   on  the  shock  are  used  in  the  interpolation  for  the  overtaken 
mesh  point.   This  is  the  only  way  the  values  on  the  shock  can 
influence  the  flow  behind  the  shock. 

Table  IV.   Percentage  error  when  a  least  squares  fit  is  used 
to  compute  n  ;s~  (    and  is  used  in  the  interpolation.   The  Lax- 
Wendroff  scheme  was  used  for  the  flow  behind  the  shock. 


Ax 

Time- 
20 

■steps 
40 

100 

0.05 
0.025 

0.98 
0.10 

0.28 
0.40 

0.66 
1.8 

Table  V.   In  this  case  the  shock  fitting  used  the  equations 
expressed  in  terms  of  p,  p  and  p.  rather  than  p,  i,      and  m, 
These  equations  are  described  in  the  report  by  Richtmyer  [3]- 
Otherwise  conditions  are  the  same  as  in  IV  above. 
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Ax 

Tlme- 

-steps 

20 

^0 

100 

0.05 

0.17 

0.27 

0.59 

0.025 

8.0 

8.5 

9.0 

0.0125 

0.12 

0.34 

0.26 

Table  VI.   In  this  case  conditions  are  the  same  as  in  V 
except  that  the  least  squares  fit  Is  not  used  for  interpolation 


Ax 


Time-steps 


20 

40 

60 

80 

100 

0.025 

0.04 

0.19 

0.06 

1.5 

0.39 

0.0125 

0.01 

0.01 

0.01 

0.01 

0.01 

Note  the  sudden  jump  in  the  error  at   80   time-steps.   It 
is  clear  that  the  least  squares  fitting  is  unsatisfactory. 

These  computations  were  carried  out  on  the  IBM  7090  at 
New  York  University.   The  author  wishes  to  thank 
Professor  Richtmyer  for  suggesting  this  problem. 
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Figure  1.   The  shock  fitting  method 
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Figure  2.   The  initial  values  of   p,   p   and  p. 
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